import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

daily_data = pd.read_csv('dataset/hn_daily.csv')

daily_data['datetime'] = pd.to_datetime(daily_data['datetime'])
daily_data.set_index('datetime', inplace=True)

daily_data = daily_data.sort_index(ascending=True)

daily_data.info()

<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 3927 entries, 2015-01-01 to 2025-10-01
Data columns (total 32 columns):
 #   Column            Non-Null Count  Dtype  
---  ------            --------------  -----  
 0   name              3927 non-null   object 
 1   tempmax           3927 non-null   float64
 2   tempmin           3927 non-null   float64
 3   temp              3927 non-null   float64
 4   feelslikemax      3927 non-null   float64
 5   feelslikemin      3927 non-null   float64
 6   feelslike         3927 non-null   float64
 7   dew               3927 non-null   float64
 8   humidity          3927 non-null   float64
 9   precip            3927 non-null   float64
 10  precipprob        3927 non-null   int64  
 11  precipcover       3927 non-null   float64
 12  preciptype        2885 non-null   object 
 13  snow              3927 non-null   int64  
 14  snowdepth         3927 non-null   int64  
 15  windgust          3927 non-null   float64
 16  windspeed         3927 non-null   float64
 17  winddir           3927 non-null   float64
 18  sealevelpressure  3927 non-null   float64
 19  cloudcover        3927 non-null   float64
 20  visibility        3927 non-null   float64
 21  solarradiation    3927 non-null   float64
 22  solarenergy       3927 non-null   float64
 23  uvindex           3927 non-null   int64  
 24  severerisk        1361 non-null   float64
 25  sunrise           3927 non-null   object 
 26  sunset            3927 non-null   object 
 27  moonphase         3927 non-null   float64
 28  conditions        3927 non-null   object 
 29  description       3927 non-null   object 
 30  icon              3927 non-null   object 
 31  stations          3927 non-null   object 
dtypes: float64(20), int64(4), object(8)
memory usage: 1012.4+ KB

daily_data.describe()

plt.figure(figsize=(12, 6))
plt.plot(daily_data.index, daily_data['temp'], color='red', alpha=0.7, linewidth=0.5)
plt.title('Daily Average Temperature (2015-2025)', fontsize=14)
plt.xlabel('Date')
plt.ylabel('Temperature (°C)')
plt.grid(True, alpha=0.3)

plt.figure(figsize=(12, 6))
weekly_temp = daily_data['temp'].rolling(window=7).mean()
plt.plot(daily_data.index, weekly_temp, color='blue', alpha=0.7, linewidth=0.9)
plt.title('Weekly Average Temperature (2015-2025)', fontsize=14)
plt.xlabel('Week')
plt.ylabel('Temperature (°C)')
plt.grid(True, alpha=0.3)

plt.figure(figsize=(12, 6))
monthly_temp = daily_data['temp'].resample('M').mean()
plt.plot(monthly_temp.index, monthly_temp.values, color='blue', linewidth=1, marker='o', markersize=3)
plt.title('Monthly Average Temperature', fontsize=14, fontweight='bold')
plt.xlabel('Month')
plt.ylabel('Temperature (°C)')
plt.grid(True, alpha=0.3)

C:\Users\nguye\AppData\Local\Temp\ipykernel_25352\1123071436.py:2: FutureWarning: 'M' is deprecated and will be removed in a future version, please use 'ME' instead.
  monthly_temp = daily_data['temp'].resample('M').mean()

# Create ribbon chart for Temperature and Dew Point using monthly averages
# Resample to monthly averages for cleaner visualization
monthly_temp_dew = daily_data[['temp', 'dew']].resample('M').mean()

plt.figure(figsize=(14, 7))

# Plot temperature and dew point lines
plt.plot(monthly_temp_dew.index, monthly_temp_dew['temp'], color='red', alpha=0.9, 
         linewidth=2, label='Monthly Avg Temperature', marker='o', markersize=4)
plt.plot(monthly_temp_dew.index, monthly_temp_dew['dew'], color='blue', alpha=0.9, 
         linewidth=2, label='Monthly Avg Dew Point', marker='o', markersize=4)

# Fill the area between the two lines
plt.fill_between(monthly_temp_dew.index, monthly_temp_dew['temp'], monthly_temp_dew['dew'], 
                 alpha=0.3, color='orange', label='Temperature-Dew Point Spread')

plt.title('Monthly Average Temperature and Dew Point Relationship (2015-2025)', 
          fontsize=16, fontweight='bold')
plt.xlabel('Date', fontsize=12)
plt.ylabel('Temperature (°C)', fontsize=12)
plt.legend(loc='upper left', fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

C:\Users\nguye\AppData\Local\Temp\ipykernel_25352\1532930895.py:3: FutureWarning: 'M' is deprecated and will be removed in a future version, please use 'ME' instead.
  monthly_temp_dew = daily_data[['temp', 'dew']].resample('M').mean()

# Split data chronologically for time series
# Use 70% for training, 15% for validation (dev), 15% for testing
train_size = int(len(daily_data) * 0.7)
dev_size = int(len(daily_data) * 0.15)

# Split the data chronologically
train_data = daily_data.iloc[:train_size]
dev_data = daily_data.iloc[train_size:train_size + dev_size]
test_data = daily_data.iloc[train_size + dev_size:]

print(f"Training set: {train_data.index[0]} to {train_data.index[-1]} ({len(train_data)} samples)")
print(f"Development set: {dev_data.index[0]} to {dev_data.index[-1]} ({len(dev_data)} samples)")
print(f"Test set: {test_data.index[0]} to {test_data.index[-1]} ({len(test_data)} samples)")

Training set: 2015-01-01 00:00:00 to 2022-07-10 00:00:00 (2748 samples)
Development set: 2022-07-11 00:00:00 to 2024-02-19 00:00:00 (589 samples)
Test set: 2024-02-20 00:00:00 to 2025-10-01 00:00:00 (590 samples)

# Calculate missing values for each column
missing_data = pd.DataFrame({
    'Column': train_data.columns,
    'Missing_Count': train_data.isnull().sum(),
    'Missing_Percentage': (train_data.isnull().sum() / len(train_data) * 100).round(2)
})

# Sort by missing percentage
missing_data = missing_data.sort_values('Missing_Percentage', ascending=False)

# Filter to show only columns with missing values
missing_data_filtered = missing_data[missing_data['Missing_Count'] > 0]

# Visualize missing data
fig, axes = plt.subplots(1, 2, figsize=(16, 6))

# Bar plot of missing percentages
if len(missing_data_filtered) > 0:
    axes[0].barh(missing_data_filtered['Column'], missing_data_filtered['Missing_Percentage'], 
                 color='coral', edgecolor='darkred')
    axes[0].set_xlabel('Missing Percentage (%)', fontsize=12)
    axes[0].set_ylabel('Features', fontsize=12)
    axes[0].set_title('Missing Data by Feature (Training Set)', fontsize=14, fontweight='bold')
    axes[0].grid(True, alpha=0.3, axis='x')
    
    # Add percentage labels
    for i, (col, pct) in enumerate(zip(missing_data_filtered['Column'], 
                                        missing_data_filtered['Missing_Percentage'])):
        axes[0].text(pct + 0.5, i, f'{pct}%', va='center', fontsize=10)
else:
    axes[0].text(0.5, 0.5, 'No Missing Values', ha='center', va='center', 
                 fontsize=14, fontweight='bold')
    axes[0].set_xlim(0, 1)
    axes[0].set_ylim(0, 1)
    axes[0].axis('off')

# Heatmap of missing data pattern (sample)
# Show a sample of rows to visualize patterns
sample_size = min(100, len(train_data))
sample_data = train_data.head(sample_size)
missing_matrix = sample_data.isnull().astype(int)

# Only show columns with missing values
if len(missing_data_filtered) > 0:
    cols_with_missing = missing_data_filtered['Column'].tolist()
    missing_matrix_filtered = missing_matrix[cols_with_missing]
    
    sns.heatmap(missing_matrix_filtered.T, cmap='RdYlGn_r', cbar=True, 
                yticklabels=cols_with_missing, xticklabels=False, ax=axes[1])
    axes[1].set_title(f'Missing Data Pattern (First {sample_size} Records)', 
                      fontsize=14, fontweight='bold')
    axes[1].set_xlabel('Record Index', fontsize=12)
    axes[1].set_ylabel('Features', fontsize=12)
else:
    axes[1].text(0.5, 0.5, 'No Missing Values\nto Display', ha='center', va='center', 
                 fontsize=14, fontweight='bold')
    axes[1].set_xlim(0, 1)
    axes[1].set_ylim(0, 1)
    axes[1].axis('off')

plt.tight_layout()
plt.show()

# Select numeric columns for distribution plots
numeric_cols = train_data.select_dtypes(include=[np.number]).columns.tolist()

# Remove tempmax and tempmin as they're collinear with temp
cols_to_exclude = ['tempmax', 'tempmin', 'feelslikemax', 'feelslikemin', 'precipcover', 'snow', 'snowdepth', 'severerisk', 'moonphase']
numeric_cols = [col for col in numeric_cols if col not in cols_to_exclude]

# Create subplots
n_cols = 4  # Number of columns in the grid
n_rows = int(np.ceil(len(numeric_cols) / n_cols))

fig, axes = plt.subplots(n_rows, n_cols, figsize=(20, n_rows * 4))
axes = axes.flatten()

# Plot distribution for each numeric column
for idx, col in enumerate(numeric_cols):
    sns.histplot(train_data[col], bins=30, kde=True, ax=axes[idx], color='steelblue')
    axes[idx].set_title(f"Distribution of {col}", fontsize=12, fontweight='bold')
    axes[idx].set_xlabel(col)
    axes[idx].set_ylabel("Frequency")
    axes[idx].grid(True, alpha=0.3)

# Hide unused subplots
for idx in range(len(numeric_cols), len(axes)):
    axes[idx].set_visible(False)

plt.tight_layout()
plt.suptitle("Distribution of All Numeric Weather Attributes in Hanoi", 
             fontsize=16, fontweight='bold', y=1.001)
plt.show()

# Select categorical columns
categorical_cols = train_data.select_dtypes(include=['object']).columns.tolist()

# Remove columns that are not useful for count plots (stations, source, name, sunrise, sunset)
# Keep description to analyze weather descriptions
cols_to_exclude = ['stations', 'name', 'sunrise', 'sunset', 'source']
categorical_cols = [col for col in categorical_cols if col not in cols_to_exclude]

print(f"Categorical columns to visualize: {categorical_cols}")

# Create subplots for categorical variables
n_cols = 2  # Number of columns in the grid
n_rows = int(np.ceil(len(categorical_cols) / n_cols))

fig, axes = plt.subplots(n_rows, n_cols, figsize=(16, n_rows * 5))
if n_rows == 1:
    axes = axes.reshape(1, -1)
axes = axes.flatten()

# Plot count plot for each categorical column
for idx, col in enumerate(categorical_cols):
    # Get value counts and filter out NaN
    value_counts = train_data[col].value_counts()
    
    # For description field, show only top 15 to keep plot readable
    if col == 'description':
        value_counts = value_counts.head(15)
    
    # Create count plot
    sns.countplot(data=train_data, y=col, order=value_counts.index, 
                  ax=axes[idx], palette='viridis')
    axes[idx].set_title(f"Count Plot of {col}", fontsize=12, fontweight='bold')
    axes[idx].set_xlabel("Count")
    axes[idx].set_ylabel(col)
    axes[idx].grid(True, alpha=0.3, axis='x')
    
    # Add value counts as text
    for i, (category, count) in enumerate(value_counts.items()):
        axes[idx].text(count, i, f' {count}', va='center', fontsize=9)

# Hide unused subplots
for idx in range(len(categorical_cols), len(axes)):
    axes[idx].set_visible(False)

plt.tight_layout()
plt.suptitle("Count Plots of Categorical Weather Variables in Hanoi", 
             fontsize=16, fontweight='bold', y=1.001)
plt.show()

Categorical columns to visualize: ['preciptype', 'conditions', 'description', 'icon']

C:\Users\nguye\AppData\Local\Temp\ipykernel_25352\1545306480.py:30: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `y` variable to `hue` and set `legend=False` for the same effect.

  sns.countplot(data=train_data, y=col, order=value_counts.index,
C:\Users\nguye\AppData\Local\Temp\ipykernel_25352\1545306480.py:30: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `y` variable to `hue` and set `legend=False` for the same effect.

  sns.countplot(data=train_data, y=col, order=value_counts.index,
C:\Users\nguye\AppData\Local\Temp\ipykernel_25352\1545306480.py:30: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `y` variable to `hue` and set `legend=False` for the same effect.

  sns.countplot(data=train_data, y=col, order=value_counts.index,
C:\Users\nguye\AppData\Local\Temp\ipykernel_25352\1545306480.py:30: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `y` variable to `hue` and set `legend=False` for the same effect.

  sns.countplot(data=train_data, y=col, order=value_counts.index,

# Statistical summary for each categorical variable
categorical_cols_all = train_data.select_dtypes(include=['object']).columns.tolist()
cols_to_exclude = ['stations', 'name', 'sunrise', 'sunset', 'description', 'source']
categorical_cols_summary = [col for col in categorical_cols_all if col not in cols_to_exclude]

print("="*80)
print("CATEGORICAL VARIABLES STATISTICAL SUMMARY")
print("="*80)

for col in categorical_cols_summary:
    print(f"\n{'='*80}")
    print(f"Variable: {col.upper()}")
    print(f"{'='*80}")
    
    # Get value counts
    value_counts = train_data[col].value_counts()
    total_count = len(train_data)
    missing_count = train_data[col].isna().sum()
    
    # Create a summary dataframe
    summary_df = pd.DataFrame({
        'Category': value_counts.index,
        'Count': value_counts.values,
        'Percentage': (value_counts.values / total_count * 100).round(2)
    })
    
    print(f"\nTotal Records: {total_count}")
    print(f"Non-Missing Values: {total_count - missing_count} ({((total_count - missing_count) / total_count * 100):.2f}%)")
    print(f"Missing Values: {missing_count} ({(missing_count / total_count * 100):.2f}%)")
    print(f"Unique Categories: {len(value_counts)}")
    print(f"\nTop 10 Categories:")
    print(summary_df.head(10).to_string(index=False))

# Special handling for description field (too many unique values)
print(f"\n{'='*80}")
print(f"Variable: DESCRIPTION (Text Field)")
print(f"{'='*80}")
description_counts = train_data['description'].value_counts()
total_count = len(train_data)
missing_count = train_data['description'].isna().sum()

print(f"\nTotal Records: {total_count}")
print(f"Non-Missing Values: {total_count - missing_count} ({((total_count - missing_count) / total_count * 100):.2f}%)")
print(f"Missing Values: {missing_count} ({(missing_count / total_count * 100):.2f}%)")
print(f"Unique Descriptions: {len(description_counts)}")
print(f"\nTop 15 Most Common Descriptions:")
description_summary = pd.DataFrame({
    'Description': description_counts.head(15).index,
    'Count': description_counts.head(15).values,
    'Percentage': (description_counts.head(15).values / total_count * 100).round(2)
})
print(description_summary.to_string(index=False))
print(f"\nNote: Due to high variability ({len(description_counts)} unique descriptions),")
print(f"this field is best suited for NLP analysis rather than direct categorical encoding.")

================================================================================
CATEGORICAL VARIABLES STATISTICAL SUMMARY
================================================================================

================================================================================
Variable: PRECIPTYPE
================================================================================

Total Records: 2748
Non-Missing Values: 2021 (73.54%)
Missing Values: 727 (26.46%)
Unique Categories: 1

Top 10 Categories:
Category  Count  Percentage
    rain   2021       73.54

================================================================================
Variable: CONDITIONS
================================================================================

Total Records: 2748
Non-Missing Values: 2748 (100.00%)
Missing Values: 0 (0.00%)
Unique Categories: 6

Top 10 Categories:
              Category  Count  Percentage
Rain, Partially cloudy   1693       61.61
      Partially cloudy    569       20.71
        Rain, Overcast    311       11.32
                 Clear    130        4.73
              Overcast     33        1.20
                  Rain     12        0.44

================================================================================
Variable: ICON
================================================================================

Total Records: 2748
Non-Missing Values: 2748 (100.00%)
Missing Values: 0 (0.00%)
Unique Categories: 4

Top 10 Categories:
         Category  Count  Percentage
             rain   2016       73.36
partly-cloudy-day    569       20.71
        clear-day    130        4.73
           cloudy     33        1.20

================================================================================
Variable: DESCRIPTION (Text Field)
================================================================================

Total Records: 2748
Non-Missing Values: 2748 (100.00%)
Missing Values: 0 (0.00%)
Unique Descriptions: 33

Top 15 Most Common Descriptions:
                                                               Description  Count  Percentage
                                         Partly cloudy throughout the day.    531       19.32
                               Partly cloudy throughout the day with rain.    419       15.25
Partly cloudy throughout the day with a chance of rain throughout the day.    274        9.97
  Partly cloudy throughout the day with rain in the morning and afternoon.    208        7.57
                 Partly cloudy throughout the day with early morning rain.    195        7.10
                     Partly cloudy throughout the day with afternoon rain.    192        6.99
                Partly cloudy throughout the day with late afternoon rain.    184        6.70
                       Partly cloudy throughout the day with morning rain.    133        4.84
                                      Clear conditions throughout the day.    129        4.69
                                Cloudy skies throughout the day with rain.    106        3.86
 Cloudy skies throughout the day with a chance of rain throughout the day.     85        3.09
                Partly cloudy throughout the day with rain clearing later.     78        2.84
   Cloudy skies throughout the day with rain in the morning and afternoon.     38        1.38
                                          Cloudy skies throughout the day.     33        1.20
                  Cloudy skies throughout the day with early morning rain.     29        1.06

Note: Due to high variability (33 unique descriptions),
this field is best suited for NLP analysis rather than direct categorical encoding.

features = ["temp", "dew", "humidity", "precip", "precipprob", "windgust", "cloudcover", "solarradiation", "solarenergy", "windspeed", "winddir", "sealevelpressure", "uvindex"]
corr_matrix = daily_data[features].corr()

plt.figure(figsize=(12, 10))
sns.heatmap(corr_matrix, annot=True, fmt=".2f", cmap="coolwarm", cbar=True)
plt.title("Correlation Matrix of Weather Features", fontsize=16)
plt.show()

# Create scatter plots for key variable relationships
fig, axes = plt.subplots(2, 3, figsize=(20, 12))

# 1. Temperature vs Humidity (colored by season)
train_data['month'] = train_data.index.month
train_data['season'] = train_data['month'].apply(lambda x: 'Winter' if x in [12, 1, 2] 
                                                   else 'Spring' if x in [3, 4, 5]
                                                   else 'Summer' if x in [6, 7, 8]
                                                   else 'Autumn')

for season, color in zip(['Winter', 'Spring', 'Summer', 'Autumn'], 
                         ['blue', 'green', 'red', 'orange']):
    season_data = train_data[train_data['season'] == season]
    axes[0, 0].scatter(season_data['temp'], season_data['humidity'], 
                      alpha=0.4, s=10, label=season, color=color)
axes[0, 0].set_xlabel('Temperature (°C)', fontsize=11)
axes[0, 0].set_ylabel('Humidity (%)', fontsize=11)
axes[0, 0].set_title('Temperature vs Humidity (by Season)', fontsize=12, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 2. Temperature vs Dew Point
axes[0, 1].scatter(train_data['temp'], train_data['dew'], alpha=0.3, s=10, color='purple')
axes[0, 1].set_xlabel('Temperature (°C)', fontsize=11)
axes[0, 1].set_ylabel('Dew Point (°C)', fontsize=11)
axes[0, 1].set_title('Temperature vs Dew Point', fontsize=12, fontweight='bold')
axes[0, 1].grid(True, alpha=0.3)

# 3. Cloud Cover vs Solar Radiation
axes[0, 2].scatter(train_data['cloudcover'], train_data['solarradiation'], 
                   alpha=0.3, s=10, color='coral')
axes[0, 2].set_xlabel('Cloud Cover (%)', fontsize=11)
axes[0, 2].set_ylabel('Solar Radiation (W/m²)', fontsize=11)
axes[0, 2].set_title('Cloud Cover vs Solar Radiation', fontsize=12, fontweight='bold')
axes[0, 2].grid(True, alpha=0.3)

# 4. Precipitation vs Humidity
axes[1, 0].scatter(train_data['precip'], train_data['humidity'], 
                   alpha=0.3, s=10, color='steelblue')
axes[1, 0].set_xlabel('Precipitation (mm)', fontsize=11)
axes[1, 0].set_ylabel('Humidity (%)', fontsize=11)
axes[1, 0].set_title('Precipitation vs Humidity', fontsize=12, fontweight='bold')
axes[1, 0].set_xlim(0, 100)  # Limit x-axis for better visibility
axes[1, 0].grid(True, alpha=0.3)

# 5. Wind Speed vs Wind Gust
axes[1, 1].scatter(train_data['windspeed'], train_data['windgust'], 
                   alpha=0.3, s=10, color='darkgreen')
axes[1, 1].set_xlabel('Wind Speed (km/h)', fontsize=11)
axes[1, 1].set_ylabel('Wind Gust (km/h)', fontsize=11)
axes[1, 1].set_title('Wind Speed vs Wind Gust', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3)

# 6. Temperature vs Feels Like (colored by humidity)
scatter = axes[1, 2].scatter(train_data['temp'], train_data['feelslike'], 
                            c=train_data['humidity'], cmap='YlOrRd', 
                            alpha=0.5, s=10)
axes[1, 2].set_xlabel('Temperature (°C)', fontsize=11)
axes[1, 2].set_ylabel('Feels Like (°C)', fontsize=11)
axes[1, 2].set_title('Temperature vs Feels Like (colored by Humidity)', 
                     fontsize=12, fontweight='bold')
axes[1, 2].grid(True, alpha=0.3)
plt.colorbar(scatter, ax=axes[1, 2], label='Humidity (%)')

plt.tight_layout()
plt.show()

C:\Users\nguye\AppData\Local\Temp\ipykernel_25352\222742786.py:5: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  train_data['month'] = train_data.index.month
C:\Users\nguye\AppData\Local\Temp\ipykernel_25352\222742786.py:6: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  train_data['season'] = train_data['month'].apply(lambda x: 'Winter' if x in [12, 1, 2]

# Box plots by season and weather conditions
fig, axes = plt.subplots(2, 2, figsize=(18, 12))

# 1. Temperature by Season
season_order = ['Winter', 'Spring', 'Summer', 'Autumn']
sns.boxplot(data=train_data, x='season', y='temp', order=season_order, 
            palette='Set2', ax=axes[0, 0])
axes[0, 0].set_xlabel('Season', fontsize=11)
axes[0, 0].set_ylabel('Temperature (°C)', fontsize=11)
axes[0, 0].set_title('Temperature Distribution by Season', fontsize=12, fontweight='bold')
axes[0, 0].grid(True, alpha=0.3, axis='y')

# 2. Humidity by Season
sns.boxplot(data=train_data, x='season', y='humidity', order=season_order, 
            palette='Set2', ax=axes[0, 1])
axes[0, 1].set_xlabel('Season', fontsize=11)
axes[0, 1].set_ylabel('Humidity (%)', fontsize=11)
axes[0, 1].set_title('Humidity Distribution by Season', fontsize=12, fontweight='bold')
axes[0, 1].grid(True, alpha=0.3, axis='y')

# 3. Precipitation by Season
sns.boxplot(data=train_data, x='season', y='precip', order=season_order, 
            palette='Set2', ax=axes[1, 0])
axes[1, 0].set_xlabel('Season', fontsize=11)
axes[1, 0].set_ylabel('Precipitation (mm)', fontsize=11)
axes[1, 0].set_title('Precipitation Distribution by Season', fontsize=12, fontweight='bold')
axes[1, 0].set_ylim(0, 100)  # Limit for better visibility
axes[1, 0].grid(True, alpha=0.3, axis='y')

# 4. Temperature by Top Weather Conditions
top_conditions = train_data['conditions'].value_counts().head(6).index
train_data_top_cond = train_data[train_data['conditions'].isin(top_conditions)]
sns.boxplot(data=train_data_top_cond, y='conditions', x='temp', 
            palette='viridis', ax=axes[1, 1])
axes[1, 1].set_ylabel('Weather Condition', fontsize=11)
axes[1, 1].set_xlabel('Temperature (°C)', fontsize=11)
axes[1, 1].set_title('Temperature by Weather Condition (Top 6)', fontsize=12, fontweight='bold')
axes[1, 1].grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.show()

C:\Users\nguye\AppData\Local\Temp\ipykernel_25352\24328514.py:6: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `x` variable to `hue` and set `legend=False` for the same effect.

  sns.boxplot(data=train_data, x='season', y='temp', order=season_order,
C:\Users\nguye\AppData\Local\Temp\ipykernel_25352\24328514.py:14: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `x` variable to `hue` and set `legend=False` for the same effect.

  sns.boxplot(data=train_data, x='season', y='humidity', order=season_order,
C:\Users\nguye\AppData\Local\Temp\ipykernel_25352\24328514.py:22: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `x` variable to `hue` and set `legend=False` for the same effect.

  sns.boxplot(data=train_data, x='season', y='precip', order=season_order,
C:\Users\nguye\AppData\Local\Temp\ipykernel_25352\24328514.py:33: FutureWarning: 

Passing `palette` without assigning `hue` is deprecated and will be removed in v0.14.0. Assign the `y` variable to `hue` and set `legend=False` for the same effect.

  sns.boxplot(data=train_data_top_cond, y='conditions', x='temp',

# Aggregate data by month to see patterns
monthly_stats = train_data.groupby('month').agg({
    'temp': 'mean',
    'humidity': 'mean',
    'precip': 'sum',
    'windspeed': 'mean',
    'cloudcover': 'mean',
    'solarradiation': 'mean'
}).reset_index()

fig, axes = plt.subplots(2, 3, figsize=(20, 10))

# 1. Average Temperature by Month
axes[0, 0].plot(monthly_stats['month'], monthly_stats['temp'], 
                marker='o', linewidth=2, markersize=8, color='red')
axes[0, 0].set_xlabel('Month', fontsize=11)
axes[0, 0].set_ylabel('Average Temperature (°C)', fontsize=11)
axes[0, 0].set_title('Average Temperature by Month', fontsize=12, fontweight='bold')
axes[0, 0].set_xticks(range(1, 13))
axes[0, 0].grid(True, alpha=0.3)

# 2. Average Humidity by Month
axes[0, 1].plot(monthly_stats['month'], monthly_stats['humidity'], 
                marker='o', linewidth=2, markersize=8, color='blue')
axes[0, 1].set_xlabel('Month', fontsize=11)
axes[0, 1].set_ylabel('Average Humidity (%)', fontsize=11)
axes[0, 1].set_title('Average Humidity by Month', fontsize=12, fontweight='bold')
axes[0, 1].set_xticks(range(1, 13))
axes[0, 1].grid(True, alpha=0.3)

# 3. Total Precipitation by Month
axes[0, 2].bar(monthly_stats['month'], monthly_stats['precip'], 
               color='steelblue', edgecolor='navy')
axes[0, 2].set_xlabel('Month', fontsize=11)
axes[0, 2].set_ylabel('Total Precipitation (mm)', fontsize=11)
axes[0, 2].set_title('Total Precipitation by Month', fontsize=12, fontweight='bold')
axes[0, 2].set_xticks(range(1, 13))
axes[0, 2].grid(True, alpha=0.3, axis='y')

# 4. Average Wind Speed by Month
axes[1, 0].plot(monthly_stats['month'], monthly_stats['windspeed'], 
                marker='o', linewidth=2, markersize=8, color='green')
axes[1, 0].set_xlabel('Month', fontsize=11)
axes[1, 0].set_ylabel('Average Wind Speed (km/h)', fontsize=11)
axes[1, 0].set_title('Average Wind Speed by Month', fontsize=12, fontweight='bold')
axes[1, 0].set_xticks(range(1, 13))
axes[1, 0].grid(True, alpha=0.3)

# 5. Average Cloud Cover by Month
axes[1, 1].plot(monthly_stats['month'], monthly_stats['cloudcover'], 
                marker='o', linewidth=2, markersize=8, color='gray')
axes[1, 1].set_xlabel('Month', fontsize=11)
axes[1, 1].set_ylabel('Average Cloud Cover (%)', fontsize=11)
axes[1, 1].set_title('Average Cloud Cover by Month', fontsize=12, fontweight='bold')
axes[1, 1].set_xticks(range(1, 13))
axes[1, 1].grid(True, alpha=0.3)

# 6. Average Solar Radiation by Month
axes[1, 2].plot(monthly_stats['month'], monthly_stats['solarradiation'], 
                marker='o', linewidth=2, markersize=8, color='orange')
axes[1, 2].set_xlabel('Month', fontsize=11)
axes[1, 2].set_ylabel('Average Solar Radiation (W/m²)', fontsize=11)
axes[1, 2].set_title('Average Solar Radiation by Month', fontsize=12, fontweight='bold')
axes[1, 2].set_xticks(range(1, 13))
axes[1, 2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Create pair plot for key variables
# Sample data for faster plotting (pairplot can be slow with large datasets)
sample_data = train_data.sample(n=min(1000, len(train_data)), random_state=42)

# Select key variables
key_vars = ['temp', 'humidity', 'dew', 'precip', 'cloudcover', 'windspeed']
sample_plot_data = sample_data[key_vars + ['season']].copy()

# Create the pair plot
pairplot = sns.pairplot(sample_plot_data, hue='season', 
                        palette={'Winter': 'blue', 'Spring': 'green', 
                                'Summer': 'red', 'Autumn': 'orange'},
                        diag_kind='kde', plot_kws={'alpha': 0.6, 's': 20},
                        height=2.5)
pairplot.figure.suptitle('Pair Plot of Key Weather Variables (Colored by Season)', 
                      y=1.01, fontsize=14, fontweight='bold')
plt.show()

# Outlier detection using IQR method
def detect_outliers_iqr(data, column):
    Q1 = data[column].quantile(0.25)
    Q3 = data[column].quantile(0.75)
    IQR = Q3 - Q1
    lower_bound = Q1 - 1.5 * IQR
    upper_bound = Q3 + 1.5 * IQR
    outliers = data[(data[column] < lower_bound) | (data[column] > upper_bound)]
    return len(outliers), lower_bound, upper_bound

# Select key numeric columns for outlier analysis
outlier_cols = ['temp', 'humidity', 'dew', 'precip', 'windspeed', 'windgust', 
                'cloudcover', 'solarradiation', 'visibility', 'sealevelpressure']

# Calculate outliers
outlier_summary = []
for col in outlier_cols:
    n_outliers, lower, upper = detect_outliers_iqr(train_data, col)
    outlier_summary.append({
        'Variable': col,
        'Count_Outliers': n_outliers,
        'Percentage': round(n_outliers / len(train_data) * 100, 2),
        'Lower_Bound': round(lower, 2),
        'Upper_Bound': round(upper, 2)
    })

outlier_df = pd.DataFrame(outlier_summary)
print("="*80)
print("OUTLIER ANALYSIS (IQR Method)")
print("="*80)
print(outlier_df.to_string(index=False))
print("\n" + "="*80)

# Visualize outliers with box plots
fig, axes = plt.subplots(3, 4, figsize=(20, 12))
axes = axes.flatten()

for idx, col in enumerate(outlier_cols):
    if idx < len(axes):
        bp = axes[idx].boxplot(train_data[col].dropna(), vert=True, patch_artist=True)
        for patch in bp['boxes']:
            patch.set_facecolor('lightblue')
        axes[idx].set_title(f'Box Plot: {col}', fontsize=11, fontweight='bold')
        axes[idx].set_ylabel(col)
        axes[idx].grid(True, alpha=0.3, axis='y')
        
        # Add outlier count annotation
        n_outliers = outlier_df[outlier_df['Variable'] == col]['Count_Outliers'].values[0]
        axes[idx].text(0.5, 0.95, f'Outliers: {n_outliers}', 
                      transform=axes[idx].transAxes, 
                      ha='center', va='top', fontsize=9,
                      bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

# Hide unused subplots
for idx in range(len(outlier_cols), len(axes)):
    axes[idx].set_visible(False)

plt.tight_layout()
plt.suptitle('Outlier Detection - Box Plots for Key Variables', 
             fontsize=14, fontweight='bold', y=1.001)
plt.show()

================================================================================
OUTLIER ANALYSIS (IQR Method)
================================================================================
        Variable  Count_Outliers  Percentage  Lower_Bound  Upper_Bound
            temp               2        0.07         8.90        40.90
        humidity              48        1.75        50.75       107.15
             dew              28        1.02         4.50        37.30
          precip             456       16.59        -3.29         5.48
       windspeed              37        1.35         4.00        30.40
        windgust              81        2.95         7.50        45.10
      cloudcover              42        1.53         4.50       126.90
  solarradiation               0        0.00       -96.80       400.20
      visibility              78        2.84         4.24        13.14
sealevelpressure               1        0.04       988.35      1033.55

================================================================================

# Comprehensive data quality check
print("="*80)
print("DATA QUALITY SUMMARY - TRAINING SET")
print("="*80)

# 1. Basic information
print(f"\n{'='*80}")
print("1. DATASET DIMENSIONS")
print(f"{'='*80}")
print(f"Total Records: {len(train_data):,}")
print(f"Total Features: {len(train_data.columns)}")
print(f"Date Range: {train_data.index.min()} to {train_data.index.max()}")
print(f"Time Span: {(train_data.index.max() - train_data.index.min()).days} days")

# 2. Duplicate records
print(f"\n{'='*80}")
print("2. DUPLICATE RECORDS")
print(f"{'='*80}")
n_duplicates = train_data.duplicated().sum()
print(f"Duplicate Rows: {n_duplicates} ({n_duplicates/len(train_data)*100:.2f}%)")

# 3. Missing values summary
print(f"\n{'='*80}")
print("3. MISSING VALUES SUMMARY")
print(f"{'='*80}")
total_missing = train_data.isnull().sum().sum()
total_cells = train_data.shape[0] * train_data.shape[1]
print(f"Total Missing Values: {total_missing:,} ({total_missing/total_cells*100:.2f}% of all cells)")
print(f"Features with Missing Data: {(train_data.isnull().sum() > 0).sum()}")
print(f"Features without Missing Data: {(train_data.isnull().sum() == 0).sum()}")

# 4. Data type distribution
print(f"\n{'='*80}")
print("4. DATA TYPE DISTRIBUTION")
print(f"{'='*80}")
dtype_counts = train_data.dtypes.value_counts()
for dtype, count in dtype_counts.items():
    print(f"{dtype}: {count} features")

# 5. Value range validation
print(f"\n{'='*80}")
print("5. VALUE RANGE VALIDATION")
print(f"{'='*80}")

validation_checks = []

# Check humidity (should be 0-100)
humidity_invalid = ((train_data['humidity'] < 0) | (train_data['humidity'] > 100)).sum()
validation_checks.append(('Humidity', '0-100%', humidity_invalid, 'Valid' if humidity_invalid == 0 else 'WARNING'))

# Check cloud cover (should be 0-100)
cloudcover_invalid = ((train_data['cloudcover'] < 0) | (train_data['cloudcover'] > 100)).sum()
validation_checks.append(('Cloud Cover', '0-100%', cloudcover_invalid, 'Valid' if cloudcover_invalid == 0 else 'WARNING'))

# Check precipitation (should be >= 0)
precip_invalid = (train_data['precip'] < 0).sum()
validation_checks.append(('Precipitation', '>= 0 mm', precip_invalid, 'Valid' if precip_invalid == 0 else 'WARNING'))

# Check wind speed (should be >= 0)
windspeed_invalid = (train_data['windspeed'] < 0).sum()
validation_checks.append(('Wind Speed', '>= 0 km/h', windspeed_invalid, 'Valid' if windspeed_invalid == 0 else 'WARNING'))

# Check temperature (reasonable range for Hanoi: -5 to 50°C)
temp_invalid = ((train_data['temp'] < -5) | (train_data['temp'] > 50)).sum()
validation_checks.append(('Temperature', '-5 to 50°C', temp_invalid, 'Valid' if temp_invalid == 0 else 'WARNING'))

for check in validation_checks:
    status_symbol = '✓' if check[3] == 'Valid' else '⚠'
    print(f"{status_symbol} {check[0]}: Expected {check[1]} | Invalid: {check[2]} | Status: {check[3]}")

# 6. Data completeness by feature group
print(f"\n{'='*80}")
print("6. DATA COMPLETENESS BY FEATURE GROUP")
print(f"{'='*80}")

feature_groups = {
    'Temperature': ['temp', 'tempmax', 'tempmin', 'feelslike'],
    'Moisture': ['humidity', 'dew', 'precip', 'precipprob'],
    'Wind': ['windspeed', 'windgust', 'winddir'],
    'Atmosphere': ['sealevelpressure', 'cloudcover', 'visibility'],
    'Solar': ['solarradiation', 'solarenergy', 'uvindex']
}

for group, features in feature_groups.items():
    available_features = [f for f in features if f in train_data.columns]
    if available_features:
        completeness = train_data[available_features].notna().all(axis=1).sum()
        pct = completeness / len(train_data) * 100
        print(f"{group}: {completeness:,}/{len(train_data):,} records complete ({pct:.2f}%)")

# 7. Summary
print(f"\n{'='*80}")
print("7. OVERALL DATA QUALITY ASSESSMENT")
print(f"{'='*80}")
quality_score = 100 - (total_missing/total_cells*100)
print(f"Data Completeness Score: {quality_score:.2f}%")
print(f"Data Quality Status: {'EXCELLENT' if quality_score > 95 else 'GOOD' if quality_score > 90 else 'FAIR'}")
print(f"\n{'✓ Dataset is ready for preprocessing and feature engineering!' if quality_score > 90 else '⚠ Consider addressing data quality issues before proceeding.'}")
print(f"\n{'='*80}")

================================================================================
DATA QUALITY SUMMARY - TRAINING SET
================================================================================

================================================================================
1. DATASET DIMENSIONS
================================================================================
Total Records: 2,748
Total Features: 34
Date Range: 2015-01-01 00:00:00 to 2022-07-10 00:00:00
Time Span: 2747 days

================================================================================
2. DUPLICATE RECORDS
================================================================================
Duplicate Rows: 0 (0.00%)

================================================================================
3. MISSING VALUES SUMMARY
================================================================================
Total Missing Values: 3,293 (3.52% of all cells)
Features with Missing Data: 2
Features without Missing Data: 32

================================================================================
4. DATA TYPE DISTRIBUTION
================================================================================
float64: 20 features
object: 9 features
int64: 4 features
int32: 1 features

================================================================================
5. VALUE RANGE VALIDATION
================================================================================
✓ Humidity: Expected 0-100% | Invalid: 0 | Status: Valid
✓ Cloud Cover: Expected 0-100% | Invalid: 0 | Status: Valid
✓ Precipitation: Expected >= 0 mm | Invalid: 0 | Status: Valid
✓ Wind Speed: Expected >= 0 km/h | Invalid: 0 | Status: Valid
✓ Temperature: Expected -5 to 50°C | Invalid: 0 | Status: Valid

================================================================================
6. DATA COMPLETENESS BY FEATURE GROUP
================================================================================
Temperature: 2,748/2,748 records complete (100.00%)
Moisture: 2,748/2,748 records complete (100.00%)
Wind: 2,748/2,748 records complete (100.00%)
Atmosphere: 2,748/2,748 records complete (100.00%)
Solar: 2,748/2,748 records complete (100.00%)

================================================================================
7. OVERALL DATA QUALITY ASSESSMENT
================================================================================
Data Completeness Score: 96.48%
Data Quality Status: EXCELLENT

✓ Dataset is ready for preprocessing and feature engineering!

================================================================================

ROLLING_WINDOWS = [7, 14, 21, 28, 56, 84]  # Days
LAG_PERIODS = {
    'temp': [1, 3, 7],
    'dew': [1, 3, 7], 
    'humidity': [1, 3, 7],
    'sealevelpressure': [1, 3, 7],
    'windspeed': [1, 3, 7],
    'precip': [1, 3, 7],
    'cloudcover': [1, 3, 7],
    'winddir_sin': [1, 3, 7],
    'winddir_cos': [1, 3, 7]
}
CATEGORICAL_FEATURES = ['icon']
TRAIN_RATIO = 0.7  # 70% train, 15% dev, 15% test
N_STEPS_AHEAD = 5  # Predict t+1, t+2, t+3, t+4, t+5
FEATURE_SELECTION_TOP_N = 30  # Per target

# Uncomment the following line to run the processing script
# %run src/daily_forecast_model/process.py

# From helper.py
PER_TARGET_PARAMS  # Loaded from best_params_per_target.json
MODELS_DIR = 'trained_models/'
DEVICE = 'cuda' if GPU available else 'cpu'

# Uncomment to run hyperparameter tuning (takes ~2-3 hours)
# %run src/daily_forecast_model/tune.py

# Uncomment to train final models using tuned hyperparameters
# %run src/daily_forecast_model/train.py

# Uncomment to run model evaluation
# %run src/daily_forecast_model/evaluate.py

import pandas as pd
import json
import os
from IPython.display import Image, display

# Load evaluation metrics
metrics_df = pd.read_csv('src/daily_forecast_model/evaluate_results/evaluation_metrics.csv')

print("="*70)
print("EVALUATION METRICS SUMMARY")
print("="*70)
print()
print(metrics_df.to_string(index=False))
print()
print("="*70)

# Calculate average metrics
avg_metrics = metrics_df[['MAE', 'RMSE', 'MAPE', 'R2']].mean()
print("\nAVERAGE ACROSS ALL TARGETS:")
print(f"  MAE:  {avg_metrics['MAE']:.4f}°C")
print(f"  RMSE: {avg_metrics['RMSE']:.4f}°C")
print(f"  MAPE: {avg_metrics['MAPE']:.2f}%")
print(f"  R²:   {avg_metrics['R2']:.4f}")
print("="*70)

======================================================================
EVALUATION METRICS SUMMARY
======================================================================

target      MAE     RMSE     MAPE       R2  n_samples
   t+1 1.167413 1.518291 4.944298 0.902168        585
   t+2 1.545550 1.981882 6.657903 0.833327        585
   t+3 1.646464 2.135913 7.137850 0.806368        585
   t+4 1.731595 2.206022 7.554534 0.792730        585
   t+5 1.747100 2.212801 7.577222 0.788970        585

======================================================================

AVERAGE ACROSS ALL TARGETS:
  MAE:  1.5676°C
  RMSE: 2.0110°C
  MAPE: 6.77%
  R²:   0.8247
======================================================================

from IPython.display import HTML

# Display individual scatter plots for each target in a grid layout
print("Individual Target Scatter Plots:")
print()


# Create HTML for displaying images in a row
html_content = '<div style="display: flex; flex-wrap: wrap; justify-content: space-around;">'

for i in range(1, 6):
    target = f't+{i}'
    plot_path = f'src/daily_forecast_model/evaluate_results/plots/scatter_{target}.png'
    
    if os.path.exists(plot_path):
        html_content += f'<div style="margin: 10px;"><img src="{plot_path}" width="350"/><p style="text-align: center; font-weight: bold;">{target}</p></div>'
    else:
        print(f"Plot not found: {plot_path}")

html_content += '</div>'
display(HTML(html_content))

Individual Target Scatter Plots:

# Display combined time series plot (all targets in multi-panel view)
combined_ts_path = 'src/daily_forecast_model/evaluate_results/plots/timeseries_all_targets.png'

if os.path.exists(combined_ts_path):
    print("\nCombined Time Series - All Targets (Multi-Panel):")
    display(Image(filename=combined_ts_path, width=1200))
else:
    print(f"Combined time series plot not found: {combined_ts_path}")

Combined Time Series - All Targets (Multi-Panel):

import matplotlib.pyplot as plt
import seaborn as sns

# Extract metrics for visualization
targets = metrics_df['target'].tolist()
mae_values = metrics_df['MAE'].tolist()
rmse_values = metrics_df['RMSE'].tolist()
r2_values = metrics_df['R2'].tolist()

# Create visualization
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# MAE by target
axes[0].bar(targets, mae_values, color='skyblue', edgecolor='navy')
axes[0].set_title('MAE by Forecast Horizon', fontsize=13, fontweight='bold')
axes[0].set_xlabel('Target', fontsize=11)
axes[0].set_ylabel('MAE (°C)', fontsize=11)
axes[0].grid(True, alpha=0.3, axis='y')
for i, v in enumerate(mae_values):
    axes[0].text(i, v + 0.02, f'{v:.3f}', ha='center', fontsize=9)

# RMSE by target
axes[1].bar(targets, rmse_values, color='lightcoral', edgecolor='darkred')
axes[1].set_title('RMSE by Forecast Horizon', fontsize=13, fontweight='bold')
axes[1].set_xlabel('Target', fontsize=11)
axes[1].set_ylabel('RMSE (°C)', fontsize=11)
axes[1].grid(True, alpha=0.3, axis='y')
for i, v in enumerate(rmse_values):
    axes[1].text(i, v + 0.02, f'{v:.3f}', ha='center', fontsize=9)

# R² by target
axes[2].bar(targets, r2_values, color='lightgreen', edgecolor='darkgreen')
axes[2].set_title('R² by Forecast Horizon', fontsize=13, fontweight='bold')
axes[2].set_xlabel('Target', fontsize=11)
axes[2].set_ylabel('R² Score', fontsize=11)
axes[2].set_ylim([0, 1])
axes[2].grid(True, alpha=0.3, axis='y')
for i, v in enumerate(r2_values):
    axes[2].text(i, v + 0.01, f'{v:.3f}', ha='center', fontsize=9)

plt.tight_layout()
plt.show()

# Print insights
print("\n" + "="*70)
print("PERFORMANCE INSIGHTS")
print("="*70)
print(f"\nBest performing target (lowest RMSE): {targets[rmse_values.index(min(rmse_values))]}")
print(f"Worst performing target (highest RMSE): {targets[rmse_values.index(max(rmse_values))]}")
print(f"\nRMSE range: {min(rmse_values):.4f}°C - {max(rmse_values):.4f}°C")
print(f"MAE range: {min(mae_values):.4f}°C - {max(mae_values):.4f}°C")
print(f"R² range: {min(r2_values):.4f} - {max(r2_values):.4f}")
print("="*70)

======================================================================
PERFORMANCE INSIGHTS
======================================================================

Best performing target (lowest RMSE): t+1
Worst performing target (highest RMSE): t+5

RMSE range: 1.5183°C - 2.2128°C
MAE range: 1.1674°C - 1.7471°C
R² range: 0.7890 - 0.9022
======================================================================

LAG_PERIODS_HOURS = {
    'temp': [1, 3, 6, 12, 24],        # 1h, 3h, 6h, 12h, 24h ago
    'humidity': [1, 6, 24],
    'sealevelpressure': [6, 24],
    'windspeed': [1, 3, 6, 24],
}

# All rolling windows use shift(1) before calculation
df["precip_roll_mean_6h"] = df["precip"].shift(1).rolling(6).mean()
df["humidity_roll_mean_24h"] = df["humidity"].shift(1).rolling(24).mean()
df["windspeed_roll_max_3h"] = df["windspeed"].shift(1).rolling(3).max()

# Hour of day (0-23)
df["hour"] = df.index.hour
df["hour_sin"] = np.sin(2 * np.pi * df["hour"] / 24)
df["hour_cos"] = np.cos(2 * np.pi * df["hour"] / 24)

# Night indicator (0-5 AM)
df["is_night"] = df["hour"].isin([0, 1, 2, 3, 4, 5]).astype(int)

cols_to_drop = [
    "tempmax", "tempmin",           # Daily aggregates
    "feelslikemax", "feelslikemin", # Daily extremes
    "feelslike",                     # Current hour aggregate
]

class PerHorizonWrapper:
    """Wrapper for per-horizon optimized models."""
    def __init__(self, models_dict):
        self.models_dict = models_dict  # {t+1h: model1, t+2h: model2, ...}
        self.horizons = sorted(models_dict.keys())
    
    def predict(self, X):
        # Route each horizon through its specialized model
        predictions = []
        for horizon in self.horizons:
            pred = self.models_dict[horizon].predict(X)
            predictions.append(pred.reshape(-1, 1))
        return np.hstack(predictions)

# Automatic device selection
DEVICE = 'cuda' if torch.cuda.is_available() else 'mps' if torch.backends.mps.is_available() else 'cpu'

# GPU-enabled models
xgb_params = {'tree_method': 'hist', 'device': DEVICE}
catboost_params = {'task_type': 'GPU' if DEVICE == 'cuda' else 'CPU'}

# Load and display hourly data statistics
import pandas as pd
import numpy as np

# Load hourly dataset
df_hourly = pd.read_csv('dataset/hn_hourly.csv')
df_hourly['datetime'] = pd.to_datetime(df_hourly['datetime'])
df_hourly = df_hourly.sort_values('datetime').set_index('datetime')

print("="*70)
print("HOURLY DATASET STATISTICS")
print("="*70)
print(f"\nDataset shape: {df_hourly.shape}")
print(f"Date range: {df_hourly.index.min()} to {df_hourly.index.max()}")
print(f"Total duration: {(df_hourly.index.max() - df_hourly.index.min()).days} days")
print(f"Total hours: {len(df_hourly):,}")

print("\n" + "="*70)
print("TEMPERATURE STATISTICS (Hourly)")
print("="*70)
print(df_hourly['temp'].describe())

print("\n" + "="*70)
print("DATA SPLITS")
print("="*70)
n = len(df_hourly)
n_train = int(n * 0.7)
n_dev = int(n * 0.15)
print(f"Training:   {n_train:,} hours ({n_train/24:.0f} days)")
print(f"Validation: {n_dev:,} hours ({n_dev/24:.0f} days)")
print(f"Test:       {n - n_train - n_dev:,} hours ({(n - n_train - n_dev)/24:.0f} days)")

======================================================================
HOURLY DATASET STATISTICS
======================================================================

Dataset shape: (94248, 23)
Date range: 2015-01-01 00:00:00 to 2025-10-01 23:00:00
Total duration: 3926 days
Total hours: 94,248

======================================================================
TEMPERATURE STATISTICS (Hourly)
======================================================================
count    94248.000000
mean        24.896258
std          5.613481
min          5.900000
25%         21.000000
50%         25.900000
75%         29.000000
max         41.900000
Name: temp, dtype: float64

======================================================================
DATA SPLITS
======================================================================
Training:   65,973 hours (2749 days)
Validation: 14,137 hours (589 days)
Test:       14,138 hours (589 days)

# Load and display actual model performance metrics from evaluation results
import matplotlib.pyplot as plt
import json

# Load actual evaluation results
with open('src/hourly_forecast_model/evaluate_results/evaluation_results.json', 'r') as f:
    eval_results = json.load(f)

# Extract RMSE values for all horizons
horizons = list(range(1, 25))
rmse_values = [eval_results['per_hour_metrics'][f't+{h}h']['metrics']['RMSE'] for h in horizons]
mae_values = [eval_results['per_hour_metrics'][f't+{h}h']['metrics']['MAE'] for h in horizons]
r2_values = [eval_results['per_hour_metrics'][f't+{h}h']['metrics']['R2'] for h in horizons]

# Get average metrics
avg_metrics = eval_results['average_metrics']

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(16, 10))

# RMSE by horizon
ax1.plot(horizons, rmse_values, marker='o', linewidth=2, markersize=6, color='#E74C3C')
ax1.axhline(y=avg_metrics['RMSE'], color='#3498DB', linestyle='--', 
            label=f"Average: {avg_metrics['RMSE']:.2f}°C", linewidth=2)
ax1.fill_between(horizons, rmse_values, alpha=0.3, color='#E74C3C')
ax1.set_xlabel('Forecast Horizon (hours ahead)', fontsize=12, fontweight='bold')
ax1.set_ylabel('RMSE (°C)', fontsize=12, fontweight='bold')
ax1.set_title('RMSE by Forecast Horizon', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3, linestyle=':', linewidth=0.8)
ax1.legend(fontsize=11, framealpha=0.9)
ax1.set_xticks(range(0, 25, 3))

# MAE by horizon
ax2.plot(horizons, mae_values, marker='s', linewidth=2, markersize=6, color='#3498DB')
ax2.axhline(y=avg_metrics['MAE'], color='#E74C3C', linestyle='--', 
            label=f"Average: {avg_metrics['MAE']:.2f}°C", linewidth=2)
ax2.fill_between(horizons, mae_values, alpha=0.3, color='#3498DB')
ax2.set_xlabel('Forecast Horizon (hours ahead)', fontsize=12, fontweight='bold')
ax2.set_ylabel('MAE (°C)', fontsize=12, fontweight='bold')
ax2.set_title('MAE by Forecast Horizon', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3, linestyle=':', linewidth=0.8)
ax2.legend(fontsize=11, framealpha=0.9)
ax2.set_xticks(range(0, 25, 3))

# R² by horizon
ax3.plot(horizons, r2_values, marker='^', linewidth=2, markersize=6, color='#2ECC71')
ax3.axhline(y=avg_metrics['R2'], color='#E74C3C', linestyle='--', 
            label=f"Average: {avg_metrics['R2']:.4f}", linewidth=2)
ax3.fill_between(horizons, r2_values, alpha=0.3, color='#2ECC71')
ax3.set_xlabel('Forecast Horizon (hours ahead)', fontsize=12, fontweight='bold')
ax3.set_ylabel('R² Score', fontsize=12, fontweight='bold')
ax3.set_title('R² Score by Forecast Horizon', fontsize=14, fontweight='bold')
ax3.grid(True, alpha=0.3, linestyle=':', linewidth=0.8)
ax3.legend(fontsize=11, framealpha=0.9)
ax3.set_xticks(range(0, 25, 3))

# Model distribution
model_counts = {'CatBoost': 17, 'LightGBM': 5, 'XGBoost': 2}
colors = ['#E67E22', '#2ECC71', '#9B59B6']
ax4.pie(model_counts.values(), labels=model_counts.keys(), autopct='%1.1f%%',
        colors=colors, startangle=90, textprops={'fontsize': 12, 'fontweight': 'bold'},
        wedgeprops={'edgecolor': 'white', 'linewidth': 2})
ax4.set_title('Best Model Distribution Across 24 Horizons', fontsize=14, fontweight='bold')

plt.tight_layout()
plt.show()

print("="*70)
print("HOURLY MODEL PERFORMANCE SUMMARY (Test Set)")
print("="*70)
print(f"\nOverall Metrics:")
print(f"  Average MAE:  {avg_metrics['MAE']:.4f}°C")
print(f"  Average RMSE: {avg_metrics['RMSE']:.4f}°C")
print(f"  Average MAPE: {avg_metrics['MAPE']:.2f}%")
print(f"  Average R²:   {avg_metrics['R2']:.4f}")

print(f"\nRMSE Range:")
print(f"  Best  (t+1h):  {min(rmse_values):.4f}°C")
print(f"  Worst (t+24h): {max(rmse_values):.4f}°C")
print(f"  Spread:        {max(rmse_values) - min(rmse_values):.4f}°C")

print(f"\nModel Selection (Per-Horizon Optimization):")
for model, count in sorted(model_counts.items(), key=lambda x: -x[1]):
    print(f"  {model}: {count}/24 horizons ({count/24*100:.1f}%)")

======================================================================
HOURLY MODEL PERFORMANCE SUMMARY (Test Set)
======================================================================

Overall Metrics:
  Average MAE:  1.1912°C
  Average RMSE: 1.6070°C
  Average MAPE: 4.85%
  Average R²:   0.9068

RMSE Range:
  Best  (t+1h):  0.7922°C
  Worst (t+24h): 1.9692°C
  Spread:        1.1770°C

Model Selection (Per-Horizon Optimization):
  CatBoost: 17/24 horizons (70.8%)
  LightGBM: 5/24 horizons (20.8%)
  XGBoost: 2/24 horizons (8.3%)

feature	meaning	values	formula
tempmax	Maximum daily temperature	°C	Max(temp)
tempmin	Minimum daily temperature	°C	Min(temp)
temp	Average daily temperature	°C	(tempmax + tempmin) / 2
feelslikemax	Maximum daily feels-like temperature	°C	Calculated using Heat Index (combination of temperature and humidity)
feelslikemin	Minimum daily feels-like temperature	°C	Calculated using Wind Chill (combination of temperature and wind speed)
feelslike	Average daily feels-like temperature	°C	Heat Index (hot and humid) or Wind Chill (cold, windy)
dew	The dew point is the temperature at which air must be cooled to become saturated with moisture and start forming dew	°C
humidity	Daily values of humidity are the mean of the hourly values (average humidity)	0–100 (%)	Daily values of humidity are the mean of the hourly values
precip	Daily precipitation amount	mm	Sum of hourly precipitation values measured in the day
precipprob	Probability of precipitation	0–100 (%)	Forecasted probability of precipitation
precipcover	Percentage of time with precipitation in the day	0–100 (%)	(Number of hours with precipitation / 24) × 100
preciptype	Type of precipitation	rain, snow, ice, freezing rain
snow	The amount of new snow that has fallen in the time period	cm
snowdepth	The average amount of snow currently on the ground for the time period	cm
windgust	The maximum wind speed measured over a short period of time (typically less than 20 seconds)	km/h or mph	Max(wind speed)
windspeed	Average wind speed	km/h	Average wind speed during the day
winddir	Wind direction	0–360° (0 degrees = North, 90 = East, 180 = South, 270 = West, back to 360 = North)
sealevelpressure	The atmospheric pressure at sea level (removes reduction in pressure due to altitude). Also called Sea Level Pressure	mb	Daily values of sea level pressure are the mean of the hourly values
cloudcover	The amount of sky that is covered by cloud (cloud cover)	0–100 (%)	Daily values include the mean of the hourly cloud coverage values
visibility	The distance that can be seen in daylight (visibility)	km	Daily values are created from the mean of the hourly values
solarradiation	Solar radiation	W/m²	Daily values of solar radiation are the mean of the hourly values
solarenergy	Solar energy during the day	MJ/m²	Daily values of solar energy are the sum of the hourly values
uvindex	UV index	0–10	The UV index is calculated based on short-wave solar radiation, cloudiness, cloud type, time of day, time of year, and altitude. Daily values represent the maximum of hourly values
severerisk	The risk of convective storms such as thunderstorms, hail, and tornadoes	0–100	Low Risk: Values less than 30 indicate a low likelihood of severe weather events; Moderate Risk: Values between 30 and 70 suggest a moderate probability of severe weather activity; High Risk: Values above 70 denote a high potential for severe weather, requiring caution and preparedness
sunrise	Sunrise time	HH:MM:SS
sunset	Sunset time	HH:MM:SS
moonphase	Moon phase	0–1	0 – new moon, 0–0.25 – waxing crescent, 0.25 – first quarter, 0.25–0.5 – waxing gibbous, 0.5 – full moon, 0.5–0.75 – waning gibbous, 0.75 – last quarter, 0.75–1 – waning crescent
conditions	Notable weather conditions reported at a particular location (general weather conditions)	thunderstorms, rainfall, etc.
description	Description of the weather for the day	Text description from API
icon	A weather icon	clear-day, clear-night, rain, snow, etc.
stations	Weather station code	VVNB
source	Data source	obs
name	Station or geographical location name of where the data was recorded	Hanoi

	tempmax	tempmin	temp	feelslikemax	feelslikemin	feelslike	dew	humidity	precip	precipprob	...	windspeed	winddir	sealevelpressure	cloudcover	visibility	solarradiation	solarenergy	uvindex	severerisk	moonphase
count	3927.000000	3927.000000	3927.000000	3927.000000	3927.000000	3927.000000	3927.000000	3927.00000	3927.000000	3927.000000	...	3927.000000	3927.000000	3927.000000	3927.000000	3927.000000	3927.000000	3927.000000	3927.00000	1361.000000	3927.000000
mean	28.515152	21.934963	24.898370	32.620550	22.703107	27.397021	20.488541	78.27344	5.226019	72.167049	...	17.394423	126.961268	1011.085358	63.828928	8.605602	155.118386	13.388719	5.84441	28.518736	0.483990
std	5.842098	4.870462	5.139137	9.530943	6.075115	7.734238	5.783574	10.76699	15.895322	44.823360	...	4.999515	97.746225	7.142712	23.636997	1.756703	77.239502	6.676342	2.47802	24.089535	0.288759
min	8.000000	5.900000	7.000000	8.000000	2.100000	4.600000	-2.000000	31.00000	0.000000	0.000000	...	6.500000	0.000000	991.100000	0.000000	2.400000	0.000000	0.000000	0.00000	5.000000	0.000000
25%	24.100000	18.000000	21.000000	24.100000	18.000000	21.000000	16.850000	72.20000	0.000000	0.000000	...	13.900000	61.900000	1005.200000	49.350000	7.700000	91.100000	7.800000	4.00000	10.000000	0.250000
50%	29.500000	23.000000	26.100000	32.700000	23.000000	27.100000	22.500000	79.80000	0.267000	100.000000	...	16.600000	109.300000	1010.700000	67.300000	9.000000	154.100000	13.300000	6.00000	10.000000	0.500000
75%	33.000000	26.000000	29.000000	41.300000	26.000000	33.900000	25.200000	85.80000	2.595000	100.000000	...	20.500000	131.900000	1016.500000	82.300000	9.800000	216.600000	18.700000	8.00000	60.000000	0.750000
max	41.900000	31.000000	35.500000	52.500000	40.400000	46.200000	28.500000	100.00000	333.117000	100.000000	...	53.600000	359.800000	1037.200000	100.000000	12.600000	323.500000	27.900000	10.00000	100.000000	0.980000

Metric	Daily Model	Hourly Model
Forecast horizons	5 days	24 hours
Average RMSE	1.88°C	1.61°C ✓
Average MAE	-	1.19°C
Best RMSE	1.52°C (t+1)	0.79°C (t+1h) ✓
Average R²	0.825	0.907 ✓
Training samples	~3,000 days	~65,000 hours
Feature count	30	30

Target	MAE (°C)	RMSE (°C)	MAPE (%)	R² Score
t+1	1.167	1.518	4.94	0.902
t+2	1.546	1.982	6.66	0.833
t+3	1.646	2.136	7.14	0.806
t+4	1.732	2.206	7.55	0.793
t+5	1.747	2.213	7.58	0.789

🌦️ Weather Features Table¶

Data Overview¶

Split the data¶

Exploratory Data Analysis¶

Univariate Analysis¶

Multivariate Analysis¶

Correlation Analysis¶

Scatter Plot Analysis¶

Box Plot Analysis by Categories¶

Temporal Patterns Analysis¶

Pair Plot for Key Variables¶

Outlier Detection¶

Data Quality Summary¶

Data Preprocessing and Feature Engineering¶

Configuration¶

Main Pipeline (Execution Order)¶

1. Load Data¶

2. Apply Feature Engineering to Full Dataset¶

3. Split Data (70/15/15)¶

4. Create Multi-Day Targets¶

5. Drop NaN Rows¶

6. Separate Features and Targets¶

7. Combined Feature Selection¶

8. Per-Target Processing Loop¶

Output Structure¶

Key Design Decisions¶

Run Preprocessing¶

Model Training and Hyperparameter Tuning¶

Two-Stage Tuning Strategy¶

Stage 1: Architecture Selection (40 trials)¶

Stage 2: Per-Target Deep Tuning (100 trials × 5 targets)¶

Final Model Training¶

Training Configuration¶

Training Process¶

Training Metadata¶

Key Design Decisions¶

ClearML Integration¶

Run Training Pipeline¶

Model Evaluation¶

Evaluation Process¶

Run Evaluation¶

Load and Display Evaluation Results¶

Scatter Plots - Actual vs Predicted¶

Time Series Plots¶

Metrics Analysis¶

Evaluation Insights¶

1. Performance by Forecast Horizon¶

2. Key Findings¶

3. Practical Implications¶

4. Model Strengths¶

5. Comparison Benchmark¶

Interface¶

Dealing with Hourly Data¶

Key Differences for Hourly Forecasting¶

1. Extended Forecast Horizon¶

2. Enhanced Lag Features¶

3. Rolling Window Statistics¶

4. Hourly Temporal Features¶

5. Data Leakage Prevention¶

6. Per-Horizon Hyperparameter Tuning¶

7. Multi-Output Model Architecture¶

8. Computational Optimizations¶

Data Statistics¶

Model Performance¶

Key Insights¶

Performance Trends¶

Model Selection Patterns¶

Comparison to Daily Model¶

Error Analysis¶

Step 7 — When to Retrain the Model¶

1. Performance Degradation Across Horizon¶

2. Data Drift¶

3. Out-of-Distribution Events¶

4. Periodic Retraining to Maintain Stability¶

5. Incremental / Mini-Update Retraining¶

6. Combined Retraining Strategy¶

7. Conclusion¶

Step 9 — ONNX for Efficient Model Deployment¶

Introduction¶

Theoretical Background¶